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ABSTRACT 


^  ' 

A^computer  simulation  model  is  developed  treating^wave 
propagation  in  a  random,  inhomogeneous  ocean  as  transmission 
through  a  linear  time-invariant,  space-variant  random  communi¬ 
cation  channel.  The  ocean  volume  is  modelled  by  an  index  of 
refraction  which  is  decomposed  into  a  depth-dependent  deter¬ 
ministic  part  and  a  depth-independent  Gaussian  zero-mean 
random  part.  Computer  simulated  output  electrical  signals 
were  generated  that  depend  on  the  complex  frequency  spectrum 
of  the  transmitted  electrical  signal,  the  far-field  beam 
pattern  of  the  transmit  array  and  the  random  transfer  function 
of  the  ocean  medium.  Output  was  generated  for  different  test 
cases.  In  all  cases  the  transmit  electrical  signal  was  repre¬ 
sented  by  a  finite  Fourier  series  and  random  cases  were 
modelled  by  a  random  number  generator.  The  computer  simulated 
output  electrical  signals  were  then  processed  by  a  3-D  DFT 
beamformer  and  the  results  for  the  deterministic  inhomogeneous 
and  random  inhomogeneous  cases  were  compared  to  the  homogeneous 
non-random  case  in  order  to  study  the  effects  of  the  medium 


on  signal  distortion  and  source  localization.  _ 
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INTRODUCTION 


I . 

Based  on  the  linearized  acoustic  wave  equation  for  small 
amplitude  acoustic  signals,  sound  propagation  through  the  ocean 
medium  can  be  modelled  as  a  signal  passing  through  a  linear 
time-variant,  space-variant  random  filter.  The  term  "space- 
variant"  implies  that  the  sound  speed  profile  is  a  function 
of  position.  The  space  variant  property  results  in  scatter 
or  angular  spread  due  to  refraction.  If  the  filter  is 
"space-variant,"  then  an  isospeed  medium  is  implied.  As  a 
result,  there  will  be  no  refraction,  and  hence,  no  scatter  or 
angular  spread  since  the  sound  rays  will  be  travelling  in 
straight  lines.  The  term  "time-variant"  implies  motion  of, 
for  example,  discrete  point  scatterers,  the  ocean  surface, 
and  the  transmit  and  receive  arrays  (apertures) .  The  time- 
variant  property  results  in  both  Doppler  spread  and  spread 
in  time  delay. 

By  using  a  linear  system's  theory  approach,  different 
propagation  paths,  such  as  direct  paths,  bottom  reflected 
paths  and  surface  reflected  paths,  can  be  modelled  as  the 
outputs  from  linear  filters.  Furthermore,  different  trans¬ 
mit  signals  and  transmit  and  receive  far-field  directivity 
functions  can  easily  be  coupled  to  various  transfer  functions 
of  the  random,  inhomogeneous  ocean  medium  in  order  to  study 
their  effects  on  signal  distortion  and  source  localization 
using  various  space-time  signal  processing  algorithms. 


Work  based  on  this  approach  has  been  done  by  Laval  [Refs. 
5,6],  Laval  and  Labasque  [Ref.  7],  Middleton  [Refs.  8,9] 
and  Ziomek  [Ref.  1] .  Middleton  [Refs.  8,9]  described  propaga¬ 
tion  phenomena  in  a  random  inhomogeneous  ocean  with  the  use 
of  space-time  operators  but  did  not  concern  himself  directly 
with  the  derivation  of  random,  time-variant,  space-variant 
transfer  functions. 

Zioniek  [Ref.  2]  derived  a  time-invariant  space-variant 
random  transfer  function  of  the  ocean  volume  based  on  the  WKB 
approximation.  Ziomek  [Ref.  3]  also  derived  an  equation  for 
the  random,  output  electrical  signal  at  each  element  in  a 
receive  planar  array  of  complex  weighted  point  sources  in 
terms  of  the  frequency  spectrum  of  the  transmitted  electrical 
signal,  the  transmit  and  receive  arrays,  and  the  random  ocean 
medium  transfer  function. 

The  purpose  of  this  thesis  is  to  implement  on  a  computer 
the  equation  for  the  output  electrical  signal  as  derived  by 
Ziomek  [Ref.  3].  This  mathematical  computer  model  simulates 
the  effects  of  wave  propagation  through  a  random  inhomogeneous 
ocean.  The  ocean  volume  is  modelled  by  the  transfer  function 
derived  in  Ziomek  [Ref.  2]  and  incorporates  an  index  of 
refraction  which  is  a  function  of  depth.  The  index  of  refrac¬ 
tion  is  decomposed  into  a  depth-dependent  deterministic  com¬ 
ponent  and  a  depth-independent  zero  mean  Gaussian  random 
component . 

Section  II  is  devoted  to  a  theoretical  description  of  the 
computer  model.  The  notation  and  geometry  of  the  computer 


model  are  explained  and  the  fundamental  equations  on  which 
the  computer  model  is  based,  are  stated.  Also,  some  additional 
constraints  are  made  to  allow  for  simpler  implementation. 

Section  III  describes  the  numerical  techniques  used  and 
the  resulting  mathematical  procedures  which  are  implemented 
in  the  program.  Furthermore,  a  brief  outline  of  the  computer 
program  is  given. 

In  Section  IV  a  number  of  test  cases  are  defined  which 
were  used  to  test  the  computer  model.  Computer  simulated 
output  electrical  signals  were  generated  for  test  cases 
incorporating  deterministic  homogeneous,  deterministic 
inhomogeneous  and  random  inhomogeneous  medium  models.  The 
output  signals  were  processed  by  a  3-D  DFT  beamformer  and 
the  results  for  the  deterministic  homogeneous  case  were  used 
as  the  baseline  case.  Both  deterministic  and  random 
inhomogeneous  cases  were  compared  to  this  baseline  case  in 
order  to  study  the  effects  of  the  medium  on  signal  distortion 

source  localization. 
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2 

Note  that  the  factor  in  equation  (2.29)  was  moved  in 

2 

front  and  replaced  by  the  factor  1/c^.  This  is  possible 
because  the  differences  in  sound  speed  at  the  different  trans¬ 
mitter  elements  are  very  small  for  all  practical  cases  and 
the  effect  on  the  amplitude  of  H(f,i,q)  can  be  neglected. 
However,  the  dependence  of  the  initial  speed  of  sound  on  the 
index  n  of  the  transmit  array  will  be  taken  into  account  in 
the  evaluation  of  the  phase  terms  as  expressed  by  the  different 
values  of  k^  in  equation  (3.1). 

Using  the  definition  of  H(f,i,q),  the  expression  for  the 
real  output  electrical  signal  at  a  receiver  element  with 
indices  (i,q)  becomes 

KMAX 

y(t,i,q)  =  2  Re  {  I  c,  H  ( kf  ,  i  ,q)  exp  (  j  27Tkf  t)  }  (3.2) 

k=l  ^  °  ° 

The  computer  program  model  first  calculates  H(f,i,q) 
as  given  in  equation  (3.1).  Then,  with  the  use  of  H(f,i,q), 
the  output  electrical  signal  is  readily  computed  according 
to  equation  (3.2). 

The  ocean  medium  transfer  function  is  possibly  random 
and  therefore  difficult  to  integrate  numerically  so  a  change 
in  the  order  of  integration  in  equation  (3.1)  is  made  to 
achieve  a  more  suitable  form  for  implementation: 


Ill .  IMPLEMENTATION  OF  THE  MODEL 

The  main  task  of  the  computer  simulation  model  is  to 
compute  and  file  the  electrical  output  signal  at  all  elements 
of  the  receive  planar  array  according  to  equation  (2.29)  ,  The 
model  was  implemented  on  an  IBM  system/370  mainframe  computer 
using  FORTRAN.  The  program  is  computation  intensive  because 
the  evaluation  of  equation  (2.29)  involves  a  double  integral 
with  nested  summations.  Therefore,  much  attention  has  been 
given  to  program  efficiency.  This  section  discusses  the 
equations  which  were  actually  implemented  and  the  overall  logic 
of  the  program. 

A.  NUMERICAL  METHODS 

The  model  computes  the  real  output  electrical  signal  at 
all  receiver  elements  as  given  by  equation  (2.29).  Note  that 
the  double  integral  with  respect  to  (w.r.t.)  direction 
cosines  u^  and  v^  behaves  like  a  frequency  response  or  trans¬ 
fer  function  H(f,i,q) ,  which  depends  on  the  position  of  the 
receiver  elements  as  indicated  by  the  indices  (i.q).  We  define 


H(f  ,i,q) 


4  £1  / 

c2  -1^ 
o 


(N-l)/2 

I  d  R. (k  ,v^,n,q) 
n  M  no  ^ 

n=- (N-1) /2 


2  2  1/2 

exp[-jk  V  AY  ]exp[-jk  (1-u  -v  )  AZ] 
•^  -^noqn  ‘^-^n  oo 


(M-l)/2 

I 

tn=-(M-l)/2 


c  exp[-jk  AX.  u  ]dv  du 
m^-^nimo  oc 


3.1 


where  is  the  fundamental  frequency  1/T^  and  the  DC  term 
is  assumed  to  be  zero.  The  complex  Fourier  series  coefficients 
Cj^  determine  the  amplitude  and  phase  relations  among  the 
different  frequency  components.  The  integral  in  equation 
(2.22)  now  reduces  to  a  summation  and  the  expression  for  the 
real  output  electrical  signal  at  a  receiver  element  with 
indices  i  and  q  becomes: 


y(t,i,q)  = 


KMAX  2  1  1 


2  Re  {  ^  c^(kf^)  /  / 

k — 1  ""I  a 

q 


(N-l)/2 


(M-l)/2 


I 


m=-(M-l)/2 


c  exp  [-jk  u  AX.  Idv  du  exp[j2'iTkf  t]}; 
m  '^  -’noim  o  o  o 


(2.29) 


generator.  These  restrictions  make  it  possible  to  evaluate 
the  transfer  function  with  the  closed  form  expressions  given 
in  equations  (2.11),  (2.16)  and  (2.19)  for  amplitude,  deter¬ 
ministic  phase  component  and  random  phase  component  of  the 
transfer  function,  respectively.  In  the  calculation  of  the 
random  phase  terms  of  the  transfer  function,  a  different 
random  number  was  drawn  for  each  combination  of  transmit 
element  depth  receive  element  depth  y^  +  qd^ , 

thereby  assuming  that  the  random  phase  terms  for  the  different 
transmission  paths  are  uncorrelated. 

3.  Frequency  Spectrum  of  the  Transmitted  Electrical 
Signa 

The  transmitted  electrical  signal  was  represented  by 
a  finite  Fourier  series  with  fundamental  period  T^  and  maximum 
frequency  Hz: 


=  KMAX  • 


(2.27) 


where  KMAX  denotes  the  total  number  of  harmonics  in  the  signal, 
This  representation  does  not  impose  a  severe  restriction  on 
the  input  signal,  since  every  finite  energy  signal  can  be 
represented  in  the  sense  of  a  least  mean-square  error  by  a 
finite  Fourier  series.  The  expression  for  the  frequency 
spectrum  becomes: 


KMAX 


X(f) 


I  [c,  6(f  -kf  )  +  c*6(f  +kf  )] 


(2.28) 
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Y-axis.  Using  Snell's  law  we  find  that  at  the  turning  point 


depth  y^: 


c(yT)  = 


(2.24) 


The  relation  between  the  angle  and  is  given  by; 


=  cos{B^) 
o  o 


(2.25) 


Combining  equations  (2.24)  and  (2.25)  into  an  inequality 
which  has  to  be  obeyed  to  avoid  a  turning  point  results  in: 


n  (y) 


c  (y) 


A  -  v^ 


(2.261 


and  we  conclude  that  obeying  the  inequality  (2.23)  also 
guarantees  the  absence  of  turning  points.  , 

2 .  Sound  Speed  Profile 

The  results  presented  in  this  thesis  are  based  upon 
a  sound  speed  profile  which  is  assumed  to  have  a  constant 
gradient  g  and  a  random  component  of  the  index  of  refraction 
nj^  which  is  not  a  function  of  depth.  Note  that  the  equations 
derived  in  Ziomek  [Refs.  2,3]  are  more  general  and  that  the 
computer  program  described  in  this  thesis  can  be  easily  modi¬ 
fied  to  handle  a  more  complex  sound  speed  profile.  For  purpose 
of  computer  simulation,  the  random  component  of  the  index  of 

refraction  was  assumed  to  be  a  zero  mean,  Gaussian  random 

2 

nxamber  with  variance  a  and  was  modelled  by  a  random  number 


receive  arrays  are  used  to  couple  an  electrical  signal  to 
the  medium  and  vice-versa. 

The  region  of  integration  over  direction  cosine  v^ 
is  limited  by  a^.  This  expresses  the  fact  that  the  transfer 
function  is  not  valid  for  direction  cosine  values  v^  approach¬ 
ing  0  (close  to  a  turning  point) .  In  most  applications  the 
limits  of  integration  for  both  integrals  can  be  further  limited 
as  shown  in  Section  III. 

B.  CONSTRAINTS  ON  THE  MODEL 

In  this  section  the  inherent  constraints  on  the  model  will 
be  investigated.  Also,  some  additional  constraints  will  be  made 
1 .  Geometrical  Constraints 

These  constraints  concern  the  transfer  function  which 
is  only  valid  under  the  two  following  restrictions: 

a.  The  assumption  made  in  the  derivation  of  the  transfer 
function,  in  order  to  validate  a  binomial  expansion  was  [Ref.  2] 

I  [n^(y)  -l]/v^]  I  <  1  (2.23) 

Thus  before  using  the  computer  model  one  should  evalu¬ 
ate  equation  (2.23)  for  the  chosen  sound  speed  profile  and 
array  locations  to  make  sure  the  inequality  is  obeyed. 

b.  Absence  of  turning  points.  The  occurrence  of  turning 
points  in  a  particular  transmitter/receiver  geometry  depends 
on  the  initial  direction  of  propagation,  described  by  the 
angle  6  between  the  initial  propagation  vector  and  the 


where 


[ln2(y^  +qd') 


AY  =  y  -  y  +  qd'  -  nd 
qn  ■'r  ^  y  y 


AX.  =  X  -  X  +  id '  -  nd 
im  r  o  X  X 


AZ  =  z  -  z 
r  o 

X(f) :  complex  frequency  spectrum  of  the  transmitted 

electrical  signal. . 


Cn  =  c(y  +ndy)  ;  s^ed  of  sound  at  transmit 
eleSent  with  index  n. 


k  =  2TTf/c (y  +  ndy)  ;  magnitude  of  initial  propagation 
vector  at  transmit  element  with  index  n. 


c^,  d^:  complex  weights  of  the  transmit  array. 


H..(f,f  ;y)  ;  Random  time- invariant  space-variant 


M 


^transfer  function  of  the  ocean  medium. 

Note  that  H  is  a  function  of  both  indices 
q  and  n,  i.e.,  it  is  a  function  of  the  y 
coordinate  of  the  source  location  y^  +  nd 
and  receiver  location  y  +qd'.  This  can^ 
also  be  expressed  by  using  the  notation 
H(k^,VQ,n,q)  expressing  H  as  a  function  of 
wave  number  k  ,  direction  cosine  v  and 
indices  n  and"q.  ° 


The  above  equations  show  how  the  medium  transfer 
function  is  used  in  describing  the  acoustic  field  at  a  receiver 
location  as  a  function  of  the  acoustic  field  at  a  transmit 
location.  The  far-field  directivity  functions  of  transmit  and 


o 


(2.18) 


MR 


=  -  —  n,  / 

V>  J 


njj(^)dC 


Substituting  the  equation  of  nj^  for  a  constant  gradient  into 
equation  (2.18)  leads  to: 

«MR  =  ^t-'r 

o  ^ 

For  a  homogeneous  medium,  equations  (2.16)  and  (2.19)  become: 


®MD 


0 


(2.20) 


MR 


o 


(2.21) 


Equations  (2.11)  through  (2.13)  for  the  amplitude  and 
phase  terms  of  the  transfer  function  show  that  the  medium  has 
the  effect  of  angle  modulating  the  transmitted  sound  field, 
also  referred  to  as  scattering. 

3 .  Output  Electrical  Signal 

The  output  electrical  signal  at  a  receive  location 
(x^+id^,y^+qd^,z^)  is  given  by  [Ref.  3]: 


00  11 

y(t,x  +id’,y  +qd',z  )  =  /  X(f)  /  /  * 

^  X.  r  Y  r  -1 


(N“l)/2  „  2  2  1/2 

rp=-  (N-l)  /2  ■*  ^ 


(M-l)/2  2  2 

y  c  exp{-j]c  u  AX.  }dv  du  exp{ j2TTft}df;  u  +v  <  1  (2.22) 

-(^^l)/2  m^^-^noamoot^"  oo- 


(-k‘/2Trf  )  /■  n^{c)nj^(C)dc 


(2.13) 


This  transfer  function  is  valid  in  the  absence  of  turning 
points  and  when  the  following  inequality  is  obeyed: 


In^(y)  -  l|/v^  <  1 


(2.14) 


In  the  case  of  a  speed  of  sound  profile  with  a  constant 
gradient  g,  the  deterministic  component  of  the  index  of 
refraction  becomes : 


(2.15) 


Substitution  into  equation  (2.12)  results  in 


-11  +  y  -yo> 


(2.16) 


As  far  as  the  random  phase  term  is  concerned  we  will 

proceed  by  assuming  that  the  random  component  of  the  index  of 

refraction  is  not  a  function  of  depth  and  can  be  represented 

2 

by  a  Gaussian,  zero  mean  random  variable  with  variance  o  : 


np(i;)  nj^  an^jj^ 


(2.17) 


where  n^^  is  N(0,a  )  and  which  is  the  normalized  random 

component  of  the  index  of  refraction,  is  N(0,1).  When  n„ 

H 

is  not  a  function  of  depth,  equation  (2.13)  becomes; 


from  the  center  of  the  transmit  array  is,  in  general,  not 


equal  to  and  is  given  by  k^: 

k„  =  2Trf/c(y^+nd  )  (2.8) 

n  -'ey 

One  can  also  use  spatial  frequencies  (in  cycles/meter)  to 
denote  directions  of  propagation: 


v^/X; 


Wq/X 


(2.9) 


2 .  Transfer  Function 

The  medium  transfer  function  is  given  by  a  function  of 
frequency  f,  spatial  frequency  f^  and  receiver  depth  y  for 
a  given  source  depth  y^  [Ref.  2] : 


(2.10) 


where  is  the  amplitude  term  given  by 

=  (2Trf  )~^/^  =  (k  V  )”^/^  (2.11) 

M  y  o  o 


and  are  the  deterministic  and  random  phase  terms, 

respectively : 


The  speed  of  sound  in  the  medium  is  assumed  to  be 
only  a  function  of  depth  and  is  given  by  c(y).  The  reference 
speed  of  sound  c^  is  the  speed  of  sound  at  the  center  of  the 
transmit  array,  i.e.,  =  c(y^).  The  index  of  refraction 

n(y)  is  a  function  of  depth  and  is  composed  of  a  determinis¬ 
tic  component  and  a  random  component: 


n(y)  =  c^/c(y)  =  np(y)  +  (2.3) 

where  n_  is  the  deterministic  component  and  n„  is  the  random 
component  of  the  index  of  refraction.  The  initial  propagation 
vector  will  be  denoted  by  k^,  =  2iTf/c^  is  the  magnitude 

of  the  initial  propagation  vector  at  depth  y  =  y^  (i.e.,  at 
the  center  of  the  transmit  array) .  The  components  of  the 
initial  propagation  vector  along  the  x,  y  and  z  axes  are 
described  in  terms  of  the  direction  cosines  u^,  v_  and  w^: 


=  k  u 
o  o 


(2.4) 


k 


y 


k  v^ 
o  o 


(2.5) 


k 


z 


k  w^ 
o  o 


(2.6) 


and 


(2.7) 


Note  that  the  magnitude  of  the  initial  propagation  vector  at 
the  transmit  elements  which  are  located  at  a  depth  nd  meters 


describes  the  wave  propagation  between  a  transmit  planar  array 
and  a  receive  planar  array  whose  centers  are  located  at 
^^o'^o'^o^  and  (x^,y^,z^),  respectively.  The  x,  y  and  z  axes 
describe  cross  range,  depth  and  range,  respectively.  The 
transmit  arra^ ,  which  is  a  parallel  to  the  XY  plane,  consists 
of  M  xN  complex  weighted  point  sources,  where  M  and  N  are  odd. 
The  sources  are  equally  spaced  in  the  X-  and  Y-directions 
with  spacings  of  d^  and  d^  meters,  respectively.  Note  that 
in  general  d^  5^  d^.  The  complex  weights  are  assumed  to  be 
separcible  and  given  by 

c^  =  a  exp  { je^}  (2.1) 

m  m  ^  m 

for  the  weighting  in  the  X-direction  and  by 

^n  “ 

for  the  weighting  in  the  Y-direction.  The  combined  weight 

factor  for  each  transmit  element  is  c  d  ,  where  m  and  n  are 

m  n 

the  indices  describing  the  transmit  element  at  position 

(x  +md  ,y  +nd  ,z  ).  The  complex  weights  are  used  for  ampli- 
O  o  ^  o 

tude  shading  and  beamsteering  of  the  transmit  beam  pattern. 

The  receive  planar  array,  which  is  also  parallel  to  the  XY 

plane,  consists  of  M'  xW  elements,  where  M'  and  N'  are  odd. 

The  sources  are  equally  spaced  in  the  X-  and  Y-directions  with 

spacings  of  d^  and  d^  meters,  respectively,  where  in  general 

d'  d'  .  The  position  of  a  receiver  element  with  indices 
X  y 

(i,q)  is  given  by  (x^+id^,y^+qd^, z^) . 
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II.  THEORY  ON  THE  COMPUTER  SIMULATION  MODEL 

A.  MODEL  DESCRIPTION 

The  mathematical  computer  model  calculates  the  output 
electrical  signal  at  each  element  of  a  receive  planar  array 
of  point  sources  in  terms  of  the  complex  frequency  spectrum 
of  the  transmitted  electrical  signal,  the  transmit  and  receive 
arrays,  and  the  random  ocean  medium  transfer  function.  The 
theory  for  the  model  is  based  on  time-variant,  space-variant 
random  filters  as  discussed  in  Ziomek  [Ref.  1].  Ziomek  [Ref. 
2]  derived  a  time-invariant  space-variant  transfer  function 
for  a  random  ocean  volume  where  the  index  of  refraction,  which 
is  composed  of  a  deterministic  and  a  random  component,  is  a 
function  of  depth.  This  transfer  function  is  based  on  the 
WKB  approximation,  which  is  a  solution  to  the  wave  equation 
when  the  index  of  refraction  is  a  function  of  depth.  Further 
more,  Ziomek  [Ref.  3]  derived  an  equation  for  the  random 
output  electrical  signals  appearing  at  each  element  in  a 
receive  planar  array  of  complex  weighted  point  sources  in 
terms  of  the  frequency  spectrum  of  the  transmitted  electrical 
signal,  the  transmit  and  receive  arrays,  and  the  random  ocean 
transfer  function.  The  computer  model  combines  and  implements 
the  equations  derived  by  Ziomek  [Refs.  2,3]  for  the  transfer 
function  and  the  output  electrical  signals. 

1 .  Geometry  of  the  Model 

Figure  1  shows  the  geometry  of  the  communication 
channel  as  simulated  by  the  computer  model .  The  model 
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=  %  / 
a 

o  q 


(N-l)/2 

I 

n=-(N-l)/2 


exp  [- jk^v^AY^^] 


1  {M-l)/2  »  7  1/7 

/  I  c  exp[-jk  u  AX.  ]exp[-jk  (1-u  -v*^)  '  AZ]c3u  dv  (3.3) 

1  fc.  1  \  /o  ni  n  o  un  n  o  o  o  o 

-1  itF=-  -1)  /2 


In  order  to  improve  efficiency,  we  concentrate  first  on 
the  inner  integral  w.r.t.  u^,  which  we  rewrite  as 


1 

J  Sj^(f,u^)exp{-k^[u^AX.  + 


(1-ul-vh  ^/^AZ]  }du 
o  o  o 


(3.4) 


where ; 


AX. 


=  x„  -  X, 


id’ 

X 


Note  that 

(M-1)/2 

S  (f,u  )  =  I  c  exp[jk  u  md  ]  (3.5) 

°  m=-(M-l)/2  n  o  X 

is  the  transmit  far-field  directivity  function  of  a  line  array 
consisting  of  equally  spaced  point-sources  as  a  function  of 
direction  cosine  u^  and  frequency  f.  Now  equation  (3.3)  reduces 
to : 


.2  1  (N-l)/2 

H(f,i,q)  =  -j  /  ■ 

'-Q  a^  n=-(N-l)/2 

/  S^(f,u^)exp{-jk^[u^AX.  +  (l-u2-v^)^/^AZ]}du^  dv^ 


(3.6) 


In  the  case  of  simple  amplitude  weighting  in  the  x-direction 
such  as  rectangular,  triangular  or  Hamming  weighting,  we  can 
obtain  closed  form  expressions  for  Sv(f,u  ).  This  improves 
program  efficiency.  For  example,  in  the  case  of  rectangular 
amplitude  weighting 


sin[k  (u  -u,  )Md  /2] 
S  (f  u  )  =  o  b  X 

'^^o  sin[k  (u  -u.  )d  /2] 

no  D  X 


(3.7) 


where  Uj^  is  the  direction  cosine  value  to  which  the  transmit 
far-field  beam  pattern  is  steered. 

Before  presenting  more  details  on  how  H(f,i,q)  is  calcu¬ 
lated,  a  brief  overview  on  numerical  and  analytical  integration 
techniques  which  were  used  is  given. 

1 .  Integration  Techniques 

a.  Method  of  Stationary  Phase 

This  method  approximates  an  integral  of  the  form 


b 

I  =  /  S(z)  exp[jg(z)]dz 

a 


(3.8) 


where  g(z)  is  a  rapidly  varying  phase  term  compared  to  the 
slowly  varying  amplitude  term  S(z).  Officer  [Ref.  10]  shows 
that  this  integral  can  be  approximated  by  using  the  fact  that 
the  main  contribution  exists  at  the  stationary  point  z^^  which 
is  defined  by 


The  limits  of  the  integral  may  be  extended  to  infinity 
because  only  the  region  around  the  stationary  point  contributes 
to  the  integral.  Assume  that  around  the  stationary  point,  z 
is  given  by 

z  =  Zgp  +  C  (3.10) 

The  approximation  of  the  integral  is  then  given  by 

I  =  S(Zgp)  /27r/|g"(Zgp)  [ exp{  i  [g  ( z^^)  ± J]  }  (3.11) 

under  the  condition: 


VI 


<< 


(3.12) 


where  the  sign  in  the  phase  term  of  equation  (3.11)  equals 


the  sign  of  the  second  derivative  g"(z^p). 


b.  Numerical  Integration  Using  Newton-Cote's  Formulas 
This  standard  technique  [Ref.  llj  is  the  method 
used  by  the  program  to  perform  a  numerical  integration  over 
a  subinterval  of  the  total  region  of  integration.  Subinter¬ 
vals  can  be  defined  by  dividing  the  total  region  of  integra¬ 
tion  into  equal  parts  or  by  a  scheme  as  described  in  the  next 
section.  Each  subinterval  is  integrated  by  evaluating  the 
Simpson’s  rule  estimate: 


where  denote  values  of  the  integrand  evaluated  at 

different  points  of  the  subinterval,  and  the  distance  between 
these  points  is  given  by  the  step  size  h.  The  original  step 
size  h  is  equal  to  half  the  size  of  the  subinterval  to  be 
integrated.  The  step  size  is  then  halved  to  obtain  a  more 
accurate  estimate  which  takes  more  points  into  account.  This 
is  repeated  until  a  tolerance  criterion  is  met.  Using  succes¬ 
sive  Simpson's  rule  estimates,  one  can  do  a  Richardson's 
extrapolation  [Ref.  12].  This  extrapolation  uses  the  last 
two  estimates  by  Simpson's  rule,  S2  and  S^^,  where  S2  uses 

half  the  step  size  of  S^.  Both  estimates  have  a  global  error 
4 

of  order  h  .  The  extrapolated  value  becomes : 

Extrapolated  value  =  ^2  +  (S2  -S^)/15  (3.14) 

6 

with  an  error  of  order  h  .  If  two  successive  extrapolated 
values  are  within  tolerance,  the  last  extrapolated  value  is 
taken  as  the  value  of  the  integral  for  that  subinterval.  This 
scheme  creates  a  semi-adaptive  integration  procedure  because 
differen.t  subintervals  can  require  different  step  sizes.  It 
could  be  further  enhanced  but  would  then  require  more  overhead 
Another  possibility  would  be  the  use  of  a  Gaussian  quadrature 
rule  which  needs  less  function  evaluations  for  the  same  order 
of  error  but  has  the  disadvantage  that  the  points  are  not 
identically  spaced  and  the  scheme  cannot  be  made  adaptive  in 
a  computation  efficient  manner. 


c.  Numerical  Integration  Using  the  Euler  Summation 
A  problem  in  evaluating  an  integral  like  equation 
(3.4)  is  convergence.  The  real  and  imaginary  parts  of  the 
integrand  of  equation  (3.4)  are  oscillatory.  Successive 
positive  and  negative  contributions  tend  to  cancel,  but  each 
contribution  itself  is  not  negligible  and  the  integral  behaves 
oscillatory  as  a  function  of  its  limits  of  integration.  There 
fore,  one  has  to  integrate  over  a  relatively  large  region  to 
obtain  convergence.  This  is  in  conflict  with  the  desired 
computational  efficiency. 

Squire  [Ref.  13;  pp.  172-173]  describes  an  inte¬ 
gration  procedure  for  oscillatory  integrands  which  splits  the 
range  of  integration  into  parts  using  the  zero  crossings  of 
the  integrand  as  points  of  separation.  The  individual  subin¬ 
tervals  can  be  evaluated  with  the  use  of  standard  methods 
(as  described  in  III.A.l.b)  and  summed  to  form  the  value  of 
the  integral.  A  technique  like  the  Euler  summation  [Ref.  13: 
pp.  172-173]  can  be  used  to  improve  convergence  of  this  sum¬ 
mation.  The  Euler  summation  obtains  from  an  original  sequence 

^1 '  ^2 '  ^3  '  •  •  •  ' 


a  second  sequence 

bi  =  a2  +  a^,  b2  =  '**'  ^k  ^  ^k  ''' 


and  a  third  sequence 


and  so  on.  It  can  be  shown  that  the  sum  of 


1  ^ 
a,  +  +  -Q  c,  +  ...  +  — m, 

214181  2"'! 

is  the  same  as  the  sum  of  the  original  sequence.  Convergence 
however  is  much  improved  which  cuts  down  on  the  number  of 
subintervals  to  be  integrated  and  improves  efficiency. 

This  method  works  well  on  integrals  where  successive 
contributions  from  subintervals  to  the  integral  are  of  decreas¬ 
ing  magnitude  and  of  alternating  sign  as  illustrated  in  Figure 
2.  To  use  this  method  it  must  be  possible  to  solve  for  the 
zero  crossings  of  the  integrand,  which  generally  requires  a 
closed  form  expression  to  solve  for  the  roots  of  the  integrand. 

2 .  Application  of  Integration  Methods 

The  methods  introduced  in  the  previous  section  are 
used  to  define  a  procedure  to  calculate  H(f,i,q).  First, 
the  integration  w.r.t.  direction  cosine  is  considered, 
which  is  given  by 


-1 


Sx(f,Uo)exp{-jkn[UoAX^ 


.+(l-u^-v^) }du^  (3.4) 

o  o  o 


Both  the  method  of  stationary  phase  and  the  numerical 
method  (Sections  III.A.l.b  and  III.A.l.c)  using  the  Euler 
summation  are  implemented.  To  apply  the  method  of  stationary 


phase  we  define  in  accordance  with  equation  (3.8) 

■'^o  (3.15) 

The  amplitude  term  S  (f,u  )  is  the  directivity  function  which 

o 

is  indeed  slowly  varying  w.r.t.  g(u^).  We  find  the  sta¬ 
tionary  point  UOSP  by  setting  g'(UOSP)  equal  to  zero: 


g'(UOSP)  =  -k  AX.  +lc  AZ(1  -u^  -  v^)”^^^u^  =  0  (3.16) 

^  n  1  n  o  o  o 

leading  to 

AX.  1  -v^  1/2 

UOSP  =  ^ — y)  (3.17) 

1  +  (AX^/AZ)^ 

For  the  values  of  the  second  and  third  derivatives  one  finds 
g"(Uo)  =  k^AZd  -u^ -v^)‘^/^(l  -  v^)  (3.18) 

g'"  (Uq)  =  3k^  AZu^d  -u^  -v^)'^/^(l  -  v^)  (3.19) 


We  now  form  the  ratio  of  the  third  and  second  derivatives 
at  u^  equal  to  the  stationary  point  UOSP: 


g'"  (UOSP) 
g*’  (UOSP) 


2  1/2^'^^?  +  AZ^) 

3AX,  d-vd'-^/^ - - - 5 - 

^  °  AZ^ 


(3.20) 


The  condition  in  equation  (3.12)  becomes 


<<  1 


(3.21) 


1  +  (AX^/AZ) 


1  -v" 


from  which  we  conclude  that  the  above  condition  is  obeyed 
only  when  AZ  >>  AX^,  which  means  that  the  method  is  only 
valid  for  relatively  small  cross  ranges  -  x^  since 
AX .  =  X  -  X  +  i  d '  . 

1  r  o  X 

The  program  implements  the  following  procedure  to 
compute  the  magnitude  and  phase  of  the  integral  w.r.t.  u^ 
[see  equation  (3.4)],  denoted  by  IWRTUO. 

1)  Calculate  the  stationary  point  UOSP 


AX 


1  -  V 


UOSP  = 


1  + (AX^/AZ) ^ 


(3.17) 


2)  Calculate  the  value  of  the  second  derivative  DRV2ND 
at  the  stationary  point; 


DRV2ND  = 


k^AZ (1  -  UOSP 


-3/2  2, 

o  o 


(3.22) 


3)  Calculate  magnitude  and  phase  of  the  integral: 


MAG  =  S^(f,UOSP)  /2Tr/DRV2ND 


(3.23) 


PHASE  =  -k  [AX.UOSP  -t-  (1  -  UOSP^  -  v^)  -t-x 

n  1  o  4 


Note  that  the  integral  is  a  function  of  v^  and  the  indices  i 
and  n  through  k^  and  AX^.  The  dependence  of  the  magnitude 
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on  the  index  n  is  negligible,  thus  DRV2ND  is  calculated  for 
all  n  using  k^.  We  will  denote  the  integral  by  IWRTUO ( ,k^, i) . 

The  numerical  method  using  the  Euler  summation  can 
also  be  applied  to  the  integration  w.r.t.  u^ .  This  numerical 
method  is  implemented  to  verify  the  result  of  the  stationary 
phase  method  and  extend  application  of  the  program  to  larger 
cross  ranges.  The  applicability  of  the  method  can  be  shown 
by  separating  the  integrand  of  equation  (3.4)  into  a  real 
and  imaginary  part 


1 

/  (S  (f,u  )cos{-k  [u  AX.  +  (1 
'  X  o  n  o  1 


2  2.1/2.„., 
u^  -  v^)  AZJ  1 
o  o 


(3.25) 


+  jS  (f,u  )sin{-k  [u  AX.  +  (1  -  u^  -  v^)  ^'^^AZ]  } )  du 
xo  nox  oo  c 


As  shown  above,  the  point  of  main  contribution  to  the 
integral  in  equation  (3.25)  is  the  stationary  point  UOSP. 

Away  from  this  point  the  integrand  becomes  more  and  more 
oscillatory  and  consequently  more  difficult  to  integrate 
numerically.  One  would  therefore  like  to  use  a  scheme  in 
which  the  size  of  the  subintervals  is  dependent  on  the  behavior 
of  the  integrand.  The  method  described  in  III.B.l.c  uses 
the  zero  crossings  of  the  integrand  as  points  of  separation, 
thus  making  sure  that  the  integrand  within  a  subinterval  is 
"well-behaved"  and  can  be  efficiently  integrated  by  standard 
numerical  methods.  This  division  into  "well-behaved"  subin¬ 
tervals  could  also  be  obtained  by  using  the  relative  extremes 
as  points  of  separation.  The  integrand  of  equation  (3.25) 


has  closed  form  solutions  to  the  zero  crossings  and  relative 
extremes  of  both  the  real  and  imaginary  parts.  The  zero 
crossings  of  the  imaginary  part  of  the  integrand  are  found 
by  solving: 

-k  AX .  u  -  k  AZ(1  -  u^  =  nn  ;  n  integer  (3.26) 

n  1  o  n  o  o  ^ 


or 


u 


-AX.  (nTT)  ±AZ[(1  -v^)k^(AX^  +AZ^)  - 

1 _ on  1 _ 

k  (AX^  +  AZ^) 
n  1 


(3.27) 


These  also  are  the  locations  of  the  relative  extremes 
of  the  cosine  term  appearing  in  the  real  part  of  the  integrand. 
In  order  to  avoid  integrating  the  real  and  imaginary  parts 
of  the  integral  separately,  the  same  subintervals  are  used 
for  both  the  real  and  imaginary  parts  of  the  integral.  Points 
of  separation  for  the  subintervals  are  the  zero  crossings  of 
the  imaginary  part  of  the  integrand,  which  coincide  with  the 
relative  extremes  for  the  real  part  of  the  integrand.  As  a 
starting  point  for  the  integration,  define  the  first  subinter¬ 
val  using  the  point  of  main  contribution,  i.e.,  the  stationary 
point  UOSP . 

Fig.  3  shows  a  typical  behavior  of  the  real  and  imagin¬ 
ary  parts  of  the  integrand  around  the  stationary  point.  The 
figure  also  indicates  the  position  of  the  subintervals. 
Integration  of  the  successive  subintervals  (from  UOSP  outward), 
will  form  alternating  terms  of  decreasing  magnitude  for  both 
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Typical  Behavior  of  Integrand  in 
Eq.  (3.25)  Around  Stationary  Point 


the  real  and  imaginary  parts  of  the  integral,  making  it 
possible  to  effectively  speed  up  the  convergence  by  using  a 
complex  Euler  summation.  The  decreasing  magnitude  is  caused 
by  the  non-linear  behavior  of  the  phase  term  and  the  decreas¬ 
ing  amplitude  of  the  integrand  which  is  given  by  the  transmit 
far-field  directivity  function  S  {f,u  ).  The  program  imple- 
ments  the  following  procedure: 

1)  Calculate  the  value  of  the  stationary  point  UOSP 
according  to  equation  (3.17).  This  value  of  u  also 
denotes  the  maximum  (negative)  value  of  the  ph§se 
term. 

2)  Calculate  NZ  the  maximum  integer  multiples  of  it 
contained  in  the  maximum  (negative)  value  of  the 
phase  term  (see  Fig.  3) . 

-k  [UOSPAX.+(l-UOSP^-v^) 

NZ  =  INTEGER(  — - - i - - -  )  (3.28) 


Note  that  NZ  denotes  the  zero  crossings  of  the  imaginary 
part  of  the  integrand  nearest  to  the  stationary  point 
UOSP. 

3)  Calculate  the  values  of  Uq  for  the  zero  crossings  of 
the  imaginary  part  of  the  integrand  denoted  by  NZ 
using  equation  (3.27). 

AX  .  (NZtt)  ±AZ  [  (l-v^)k^  (AX^  +  AZ^)  -  (NZtt) 

1  _  on  1 _ _ _ 

k  (AX^  +AZ^) 
n  1 

(3.29) 

Integrate  the  subinterval  between  these  two  values 
of  UOZ. 

4)  Decrement  NZ  and  calculate  the  values  of  U0Z;j^^2' 
denoting  the  next  zero  crossings  of  the  imaginary 
part  of  the  integrand,  by  using  the  formula  given 
in  3)  above. 


5)  Define  two  subintervals,  one  on  each  side  of  the 
stationary  point,  from  zero  crossing  to  zero  crossing 
for  the  imaginary  part  and  from  extreme  to  extreme 
for  the  real  part  of  the  integrand. 

6)  Integrate  both  sub intervals  and  add  the  contribution 
to  the  total  value  of  the  integral  using  a  complex 
Euler  summation. 

7)  Check  for  convergence  by  comparing  the  difference 
between  the  last  two  values  of  the  Euler  summation 
with  a  tolerance.  If  not  convergent,  go  back  to  4) 
and  repeat  the  procedure. 

One  more  note  can  be  made  concerning  the  numerically 
calculated  value  of  the  integral  w.r.t.  u  for  different  k  . 
Recalculation  of  the  whole  integral  for  different  values  of 
n  would  be  very  time  consuming.  We  therefore  introduce  a 
correction  term  to  obtain  values  of  the  integral  for  n  ^  0 , 
based  on  the  value  of  integral  for  n  =  0.  This  correction 
term  is  only  of  importance  to  the  phase  of  the  integral,  the 
effect  of  different  values  for  k^  on  the  magnitude  is 
neglected.  For  n  =  0  the  stationary  phase  method  results  in 
the  following  expression  for  the  phase: 


PHASE  (n  =  0)  =  -k^  [AX^U0SP+ (I-UOSP^-Vq)  ^  (3.30) 


Brekhovskikh  [Ref.  14]  demonstrates  that  taking  into 
account  higher-order  terms  in  the  method  of  stationary  phase 
does  change  the  amplitude  but  not  the  phase  of  the  resulting 
approximation.  This  shows  that  the  point  of  main  contribu¬ 
tion,  i.e.,  the  stationary  point  determines  the  phase  of  the 
integral.  If  the  method  of  stationary  phase  is  no  longer 


39 


IV.  COMPUTER  SIMULATED  DATA 

This  section  describes  the  validation  of  the  model.  It 
presents  and  interprets  the  data  generated  by  the  computer 
model  for  a  number  of  test  cases. 

A.  PRESENTATION  OF  COMPUTER  SIMULATED  DATA 

Data  generated  by  the  computer  model  consists  of  the 
frequency  response  matrix  H(f,i,q)  and  the  output  electrical 
signal  at  each  element  y(t,i,q).  Interpretation  of  the  data 
is  done  directly  by  examining  the  frequency  response  H(f,i,q) 
and  by  the  use  of  a  three-dimensional  DFT  program  to  calculate 
the  Fourier  transforms  of  the  computer  simulated  output  elec¬ 
trical  signals  w.r.t.  time  and  space.  The  program  is  based 
upon  3-dimensional  DFT  beamforming  for  planar  arrays  as  des¬ 
cribed  by  Ziomek  [Ref.  4].  First  a  brief  overview  of  the 
output  from  the  3-D  DFT  program  is  given. 

The  program  calculates  the  3-dimensional  DFT  of  the 
electrical  output  signal  from  all  receiver  elements  and  is 
given  by 

M'  N'  L' 

Y  (q,r,s)  =  I  I  I  c  d  y(£,m,n)exp[-j27TqVL]  • 

^  m=-M'  n=-N'  £=-L'  ^ 

exp  [  j  2TTrm/  (M-t  Mx)  1  exp  [  j  2iTsn/  (N+NZ)  ]  (4.1) 
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DO  100  J  =  1,  N 


100  COEF(NEW, J+1)  =  COEF(NEW,J)  +  COEF(OLD,J) 

FACTOR  =  FACTOR  /  2 

EULSUM  =  EULSUM  +  COEF (NEW , N) /FACTOR 
C  Prepare  for  next  contribution 

SWAP(  NEW,  OLD  ) 

The  program  implements  this  algoritm  in  complex  form  to 
perform  N'  summations  simultaneously. 

More  details  on  the  actual  programming  of  the  discussed 
numerical  methods  and  flowcharts,  data  structures  used,  etc 
is  considered  applied  programming  and  will  not  be  discussed 


in  this  thesis. 


is  described  in  III. A. 2  for  the  integration  w.r.t.  and  is 

similar  to  the  procedure  used  for  the  numerical  integration 

w.r.t.  u  in  Section  III. A. 2. 
o 

Implementation  of  the  integration  over  a  subinterval  using 
the  Newton  Cote's  formulas  is  straightforward  and  will  not 
be  examined  in  detail.  The  framework  for  this  part  of  the 
program  can  be  found  in  Squire  [Ref.  13:  p.  281]. 

The  Euler  summation  is  used  to  speed  up  convergence  of 
the  numerical  integration.  The  objective  is  to  form  the  sum 


from  the  original  sequence  a^^ ,  a2 /a^ » ,  •••>  a^  which  are  the 
contributions  from  the  subintervals  1  through  n,  where 

\  ^  ^k+1'  ^k  ^  ^k+1'  - 

The  algorithm  used  is  described  below  and  uses  an  array 

COEF  to  store  the  coefficients  a,b  wC  o,d  -,  ...  etc. 

n  n-1  n-2  n-3 

The  size  of  this  array  determines  how  many  contributions  can 
be  handled. 

DIMENSION  C0EF(2,SIZE) 

REAL  EULSUM, FACTOR 
INTEGER  NEW,  OLD 
C  Add  contribution  N 

COEF  (NEW,1)  =  N  th  CONTRIBUTION 


DQ  \ 

(for'  all  frequency 
components  > 


Ret  urn 


integrations  and  keep  track  of  N'  Euler  suininations  simul¬ 
taneously.  It  does,  however,  also  reduce  some  overhead. 
Convergence  checking  on  the  integration  w.r.t.  v^  is  done 
only  on  the  two  outermost  elements  with  indices  q  =  (N'-l)/2 
and  q  =  -(N'-l)/2. 

The  calculation  of  those  terms  that  multiply  IWRTUO(k  ,v  ,i) 
in  the  summation  given  by  equation  (3.46)  is  simplified  by 
precomputing  part  of  the  phase  components  of  the  ocean  medium 
transfer  function  given  by  equations  (3.36)  and  (3.37).  The 
program  does  this  in  the  subroutine  HMVKB  and  stores  the 
resulting  phase  integral  values  in  an  array.  To  obtain  a 
particular  value  for  the  phase,  the  program  only  needs  to 
multiply  with  the  factor  Especially  in  the  case  of  a 

more  complicated  sound  speed  profile,  the  computation  of 
the  phase  integrals  can  be  involved  and  precomputation  im¬ 
proves  efficiency.  The  stated  time-invariance  of  the  model 
also  forces  precomputation  since  the  random  numbers  used  to 
model  the  phase  term  9^^^  can  be  drawn  only  once.  The  random 
numbers  were  obtained  from  the  IMSL  routine  GNNML  which  was 
used  to  generate  a  series  of  N  xN'  random  numbers  with  distri¬ 
bution  N(0,1).  From  these  random  numbers  the  different  random 
phase  terms  for  all  combinations  (n,q)  are  readily  com¬ 

puted  using  equation  (3.37).  Figure  5  shows  a  simple  flow¬ 
chart  of  the  CALCHF  module. 

The  innermost  block  which  actually  calculates  H(f,i,q) 
is  shown  in  more  detail  in  Fig.  6.  The  diagram  shows  how  the 
subintervals  are  defined  by  the  program.  The  procedure  followed 


Most  important  numerical  calculations  required  by  the 
model  are  done  within  the  module  CALCHF.  The  module  CALCHF 


calculates  the  frequency  response  matrix  H(f,i,q),  i.e.,  the 
frequency  response  for  all  receiver  elements  (i,q).  This 
mainly  involves  evaluation  of  the  do;ible  integral  w.r.t.  u^ 
and  v^.  It  has  two  versions  dependent  on  the  integration 
technique  used  to  evaluate  the  integral  w.r.t.  u^ : 

CALCHF_1  Use  stationary  phase  method  to  integrate 

w.r.t.  u  . 

o 

CALCHF_2  Use  numerical  technique  to  integrate  w.r.t.  u^. 

To  avoid  recalculation  and  improve  program  efficiency 
the  outer  integration  w.r.t.  v^  is  done  simultaneously  for 
all  elements  with  the  same  index  i.  This  is  particularly 
important  when  using  the  more  time  consuming  numerical  method 
to  integrate  w.r.t.  u^.  This  can  be  done  by  separating  the 
summation  inside  the  integrand  into  a  part  dependent  on  q 
and  a  part  dependent  on  i.  Recall  that  equation  (3.35) 
gives  the  form  of  the  integrand  as  follows: 


(N-l)/2 

y  cl  H..(k  ,v  ,n,q)exp{-jk  v  AY  }lWRrU0(v  ,k  ,i) 
n=-(N-l)/2^^  no''^  ^  n  o  qn  on 


(3.46) 


So  given  a  value  for  IWRTUO  for  a  given  index  i  we  can  evalu¬ 
ate  the  integrand  for  all  values  of  q,  thus  avoiding  recalcu¬ 
lating  the  most  computation  intensive  part.  This  scheme  does 
result  in  more  complex  coding.  The  program  has  to  do  N' 


Fig.  4  Flowchart  of  Main  Program  Module 
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B .  PROGRAM  OUTLINE 


The  flowchart  in  Fig.  4  gives  an  overview  of  the  program 
logic.  This  logic  is  controlled  by  flags  which  are  read 
from  the  input  data  file  together  with  all  program  parameters. 
The  program  parameters  determine: 

-  Model  geometry:  x  ,  y  ,  7.  ,  x  ,  y  ,  z  ,  M,  N,  M' ,  N', 

d^,  d  ,  d' ,  d' .  °  °  °  ^  ^  ^ 

X  y'  X  y 

Speed  of  sound  profiles:  c^,  g,  o. 

Beam  steering  angles  for  the  transmit  array  and  the 
frequency  to  be  used  for  the  beam  steering  calculations. 

The  fundamental  period  Tq  of  the  electrical  input  signal 
and  the  total  number  of  samples  L  in  the  electrical 
output  signal.  The  program  determines  the  sampling 
period  Tg  =  Tq/L  and  calculates  the  output  electrical 
signal  at  time  instants 


-(L-l)Ty2,  0,  ...,  (L-l)T^/2 


The  values  of  the  frequency  components  in  the  electrical 
input  signal  and  their  complex  Fourier  series  coefficients 

^k- 

Tolerance  values  ERRORU  and  ERRORV  for  the  numerical 
integrations  w.r.t.  u  and  v  .  The  user  should  validate 
his  choice  for  the  toierance°values  by  comparing  them 
against  the  absolute  values  calculated  by  the  program 
for  the  integrations. 

Processing  is  done  in  accordance  with  flag  settings  and 

allows  for  a  number  of  options: 

The  frequency  response  H(f,i,q)  can  be  calculated  or 
read  from  a  previously  created  file. 

The  calculated  frequency  response  H(f,i,q)  can  be 
filed  for  later  use  and/or  printed. 

The  output  electrical  signal  can  be  calculated  and 
filed  for  further  processing  by  other  prograuns . 


described  for  the  integration  w.r.t.  u^.  The  point  of  main 
contribution  is  approximately  where  the  phase  term  a (v^)  is 
stationary.  This  point  is  used  to  begin  the  numerical  inte¬ 
gration.  The  stationary  point  VOSP  is  given  by 

VOSP  =  AY/(R^.  +  (3.43) 


which  is  a  solution  of  a' (VOSP)  =  0  where  R  .  is  the  radial 

ri 

distance  from  the  center  of  the  transmit  array  to  a  receiver 
element  with  index  i  and  is  given  by 

2  2 

Rri  =  (AX^  +  AZ^)  (3.44) 


The  points  of  separation  for  the  integral  can  be  found  by 


VOZ 


1,2 


-AY(mT)±R^^[k‘(AY^+Rj^)  -(mr)^] 


1/2 


(3.45) 


Another  method  using  a  brute  force  technique  was  used 
to  validate  the  above  described  procedure.  This  method 
integrates  w.r.t.  direction  cosine  v^  between  user  specified 
limits  of  integration.  It  divides  the  range  of  integration 
into  subintervals  of  equal  size  which  are  integrated  by  a 
standard  routine  as  described  above  in  III.A.l.b.  Both 
methods  agreed  numerically  although  the  method  using  the 
unequally  spaced  subintervals  and  the  Euler  summation  was 
much  more  efficient  in  terms  of  computation  time. 
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where  Ay  =  y  -y  .  If  we  let 
^  r  ^  o 


(N-l)/2 


exp{-j (k^-k^)  [iYv^+AX^03SEM-(l-UDSP^-v^)^^^AZ]}  -  llWI?nX)(v^,k  ,i)  |  (3. 40! 


and 


a{v  )  =  -k  [AYv  +UOSPAX.  +  ( 1-UOSP^-v^ )  ^-^^AZ] 

O  O  O  1  o 


a.{v)  =  -k  [AYv^  +  (1-v^)^/^  Ax^  +AZ^] 


o  o 


Then  the  integral  given  by  equation  (3.39)  becomes 


(3.41) 


1 

/  X(v^)exp[ ja(v^) ]  dv^ 

3 

q 

where  X(v^)  is  slowly  varying  compared  to  exp{ja(v^)}.  The 
phase  term  u(v^)  can  now  be  used  to  define  subintegrals  which 
are  "well  behaved."  This  is  accomplished  by  taking  the  roots 
of  sin(a(v^))  or  cos(a(v^))  as  the  points  of  separation  between 
the  different  subintervals.  Contributions  from  the  subintegrals 
are  in  general  of  decreasing  magnitude  and  alternating  sign, 
so  the  Euler  summation  can  be  used  to  effectively  speed  up 
convergence . 

The  implemented  procedure  uses  the  roots  of  sin(a(v^)) 
as  the  points  of  separation  and  is  similar  to  the  one  already 


v^.  The  expressions  within  the  square  brackets  for  and 

are  only  dependent  on  the  geometry  and  can  be  precomputed. 
A  random  number  generator  was  used  to  generate  a  series  of 
N(0,1)  distributed  numbers  to  be  used  for  n^^^^  in  equation 
(3.37)  . 

Because  of  the  complexity  of  the  total  phase  term  in 
the  integrand  of  equation  (3.35),  a  closed  form  approximation 
could  not  be  found.  The  numerical  method  used  is  also  based 
on  the  procedure  given  in  III.A.l.c.  The  integral  is  divided 
into  subintegrals  in  such  a  way  that  the  subintegrals  are 
all  "well  behaved."  The  problem  is  to  find  a  general  form 
which  approaches  the  phase  of  the  integrand.  This  form  can 
then  be  used  to  calculate  the  points  of  separation  in  order 
to  define  the  subintervals.  We  will  proceed  by  separating 
the  phase  term  of  the  integrand  into  a  major  part  and  a  minor 
slowly  varying  part.  Combining  the  phase  term  of  the  integral 
w.r.t.  u^  given  by  equation  (3.32)  with  equation  (3.35) 
results  in  the  following  expression  for  the  integration  w.r.t. 

1  (N-l)/2 

a  n=-(N-i)/2  ^ 

exp{-j(k  -k  )  [AYv  +AX.UOSP+(l-UOGP^-v5^'^^AZ]  }  • 
n  o  o  1  o 


lB^RrU0(v  ,k  ,i)  lexp{-jk  [v  AY+UOSPAX.  +  (l-UOSP^-y?) AZ] }dv_ 
'  on'  ^^oo  1  o  o 


(3.39) 


where  we  used  =  1500  m/sec  and  g  =  0.017  1/sec.  An  upper 

limit  for  the  correction  term  is  (k  -k  ) AZ  if  AZ  >>  AX. 

n  o 

Having  established  two  methods  to  implement  the  calcu¬ 
lation  of  the  inner  integral  w.r.t.  direction  cosine  u^,  the 
outer  integral  w.r.t.  direction  cosine  is  now  considered. 

The  integral  can  be  written  as  [see  equation  (3.6)]: 

1  (N-l)/2 

a  n=-(N-l)/2  nmno  noqn  on 

(3.35) 

The  ocean  medium  transfer  function  values  are  calculated  with 
the  use  of  equations  (2.11),  (2.16)  and  (2.18).  Note  that  the 
deterministic  and  random  phase  terms  pertain  to  the  case  of 
a  sound  speed  profile  with  a  constant  gradient.  In  applying 
these  formulae,  we  have  to  take  into  account  that  the  source 
depth  for  a  transmit  element  is  not  but  yj^+ndy.  This  results 
in  the  following  formulae  which  were  actually  implemented: 


A. .  =  ( k  V  ) 

M  o  o 


-1/2 


o  ^  ^r  ^  y 


(3.36) 

(3.37) 

(3.38) 


We  note  that  both  the  amplitude  and  phase  terms 
depend  on  the  value  of  the  wavenumber  k  and  direction  cosine 


valid  we  can  write  a  general  form  for  the  phase  based  on 
equation  (3.30) : 


PHASE  (n=0)  =  -k  [AX .  UOSP+ ( 1-UOSP^-v^)  +  t  +  A(j) 

o  1  o  4  ^ 

(3.31) 

where  A(})  represents  some  unknown  correction  term.  This  relies 
on  the  fact  that  the  stationary  point  UOSP  is  always  very 
close  to  the  point  of  main  contribution.  Based  on  this,  the 
phase  term  for  n  0  can  be  written  as: 

PHASE(n)  =  -[k  +(k  -k  )  ]  [AX.UOSP+ (1-UOSP^-v^)  ^/^AZ]+x  +  A(t)  (3.32) 
o  n  o  1  o  4 

where 


(k^-k^) [AX.UOSP  +  (1-UOSP^-v^) ^/^AZ]  (3.33) 

n  o  1  o 

represents  the  correction  term  which  is  readily  computed  for 
different  values  of  n.  This  phase  correction  term  has  to 
be  added  to  the  phase  for  n  =  0  to  obtain  phase  values  for 
n  fi  0 .  Note  that  this  correction  term  is  only  important  for 
larger  ranges  AZ,  assuming  that  AZ  >>  AX.  This  can  be  seen 
by  calculating  a  typical  value  for  where  the  element 

spacing  is  taken  as  A/2  and  the  sound  speed  profile  has  a 
constant  gradient  g: 

k  -k  =  2Trf(-^ - -)  Z  -  n  =  3.6  x  10~^  n  (3.34) 

no  c  c  c 


where : 


L;  total  number  of  time  samples  in  one  fundamental 
period  of  the  signal. 

L'  =  (L-l)/2 

:  sampling  period 

M:  total  odd  number  of  elements  in  the  x-direction 

in  the  receive  array 

M'  =  (M-l)/2 

N:  total  odd  niomber  of  elements  in  the  y-direction 

in  the  receive  array 

N'  =  (N-l)/2 

spacing  in  x-  and  y-directions ,  respectively 

complex  separable  weights  for  the  array 
given  by 

=  a  exp{j0^} 
m  m 

=  b^expi 

q,r,s:  binnumbers  of  the  DFT  which  are  related  to 
values  of  frequency  f,  direction  cosine  u 
and  direction  cosi.-j  v,  respectively. 

y(Jl,m,n);  output  electrical  signal  at  time  instant 

2.T  at  element  (m,n)  of  the  receive  array, 
s 

MZ ,  NZ :  Number  of  zeros  used  for  padding  to  increase 
the  resolution  of  the  direction  cosines  u 
and  V,  respectively. 


Note  that  the  symbols  used  to  characterize  this  planar  receive 
array  are  not  the  same  symbols  used  to  describe  the  receive 
array  in  the  computer  model  {II. A. 1).  Rectangular  amplitude 
weighting  without  beam  steering  is  used  for  all  the  results 
presented  in  this  chapter. 


In  II. B. 3  we  made  the  constraint  to  represent  the  transmit 
electrical  signal  by  a  finite  Fourier  series  with  maximum 
frequency  KMAX*f^.  Thus,  to  avoid  aliasing,  we  must  sample  the 
received  signal  with  a  sample  rate  equal  to  or  greater  than 
the  Nyquist  rate: 

f^  >  2  KMAX  f  ;  T  >  ^  KMAX  T  (4.2) 

s  —  o  o  —  s 

The  total  number  of  samples  L  to  be  taken  in  order  to  avoid 
aliasing  must  obey  the  following  inequality: 

L  ^  2  KMAX  +1  (4.3) 

The  sampling  period  is  determined  by: 

Tg  =  T^/L  (4.4) 

The  program  evaluates  the  3-dimensional  summation  in 
three  steps.  The  first  step  is  to  perform  a  DFT  w.r.t.  time 
at  each  receiver  element  (m,n) .  If  the  sound  field  incident 
upon  the  array  is  a  general  plane  wave,  then  the  received 
electrical  signal  can  be  expressed  by 

y  (  «,  ,m,n)  =  g(£T  -  t  ) 

so 

where  t  =  (u  md  +  v  nd  ) /c  denotes  the  time  delay  between 
o  o  X  o  y 

the  signals  at  the  different  elements.  The  DFT  w.r.t.  time 


becomes 


Yg  (q,iii,n) 


gUVt^JWj 


(4.5) 


c  d 
m  n 


Jl=-L' 


where  is  the  factor  exp{- j 2iTq£/L} .  Using  the  above  stated 

Li 

fact  that  the  received  signal  is  given  by  a  finite  Fourier 
series,  we  obtain 


Yg{q,in,n)  =  Cj^d^Lc^exp{- j  2-TrqfQ  (u^md^  +  v^nd^) /c }  (4.6) 

A  normalized  expression  is  defined  by 
=  Yg  (q,m,n, ) /(Cj^d^L) 

=  c  exp{-j2iTqf  (u  md  +v  nd  )/c}  (4.7) 

>5  o  o  X  o  y 

So  Yg^(q,m,n)  allows  us  to  estimate  the  amplitude  of  the  fre¬ 
quency  components  in  the  signal  at  each  element  (m,n) . 

The  next  DFT  calculation  is; 

m ' 

Y„(q,r,n)  =  I  Y  (q ,m,n) exp{ j 2umr/ (m+mZ ) }  (4.8) 

m=-M '  ^ 

where  we  have  padded  with  NZ  zeros  to  increase  resolution. 
Substitution  equation  (4.6)  into  equation  (4.8)  leads  to 


Yg (q,r,n) 


=  L  c 


qdj,  exp{-j2Tr  qf^  v^n  dy/c}  • 

c  exp{-j27Tqf  u  md  /c}w5^^,,„ 
m  ^  J  ^  o  o  X  M+MZ 


(4.9) 


and  a  normalized  expression  can  be  written  as 

YsN^'5.r,n)  =  Cg  exp{-j  2tt  q  n  d^/c}  • 

M'  M' 

I  c  exp{- j2Tr  q  f  u  m  d  /c}w5^„„/  7  a  (4.10) 

ni=-M'  ^  ^  o  o  x^  m 

This  expression  is  directly  proportional  to  the  Fourier 
series  coefficients  c  and  to  the  normalized  far-field  beam 

q 

pattern  at  a  frequency  q'f^  of  a  line  array  consisting  of  M 
elements  which  is  lying  along  the  X-axis.  The  planar  array 
has  N  such  line  arrays  as  indicated  by  the  index  n.  If  no 
beam  steering  is  done  (c  =  a  ),  then  Y(,^(q,r,n)  has  a  maximum 
approaching  c  at  the  bin  number  related  to  the  direction 

q 

cosine  value  which  is  closest  to  u  .  The  DFT  bin  number  r 

o 

is  related  to  the  direction  cosine  value  by 

u  =  rX/[ (M+MZ)d^]  ;  A  =  c/qf  (4.11) 

This  value  of  u  is  an  estimate  of  u^.  Therefore,  for  each 

frequency  component  of  a  plane  wave  propagating  in  the  u^,  v^ 

direction,  the  expression  given  by  equation  (4.10)  can  be  used 

to  estimate  the  magnitude  of  the  frequency  component  and  the 

direction  cosine  u  .  An  estimate  of  v_can  be  obtained  in  the 

o  o 

same  way  by  evaluating 


^  /c}W=;„2/  J  b  (4.12) 

n=-N  ^  n=-N 

The  above  mentioned  estimates  of  c  ,  u  and  v  are  provided 

q  o  o 

by  the  program  in  the  form  of  printouts  of  Y  (q,m,n) , 
Ygj^(q,r,n)  and  Ygj^(q,m,s).  The  printouts  are  done  for  user 
specified  values  of  m  and  n.  No  padding  with  zeros  is  done 
for  the  transform  w.r.t.  time.  Padding  for  the  spatial  trans¬ 
forms  is  done  to  obtain  a  user  specified  step  size  in  direction 
cosine  u  or  v.  Graphical  output  for  estimation  of  u^  and  v^ 
was  also  obtained.  The  complete  three-dimensional  transform 
Yg(q,r,s)  was  used  to  produce  a  three-dimensional  plot  to  aid 
in  interpretation  and  for  illustrative  purposes.  The  esti¬ 
mates  of  u^  and  v  can  also  be  used  to  calculate  estimates 
o  o 

of  the  spherical  angles  denoting  the  direction  of  the  sound 
source  as  measured  from  the  center  of  the  receive  array. 

These  angles  are  also  referred  to  as  target  bearing  and  eleva¬ 
tion  angles.  The  estimates  are  given  by 

Ip  =  tan  ^(v  /u  )  (4.13) 

o  o  o 

=  sin"^(  [u^  +  (4.14) 

where  p  and  9  are  estimates  of  the  target  bearing  and  ele- 
o  o 

vation  angles,  respectively.  Note  that  an  error  in  the 


estimation  of  either  Uq  or  Vq  or  both  affects  both  the  target 
bearing  angle  and  target  elevation  angle  estimates. 


B .  TEST  CASES 

Most  test  cases  run  by  the  model  are  kept  simple  to  aid  in 
interpretation.  All  cases  are  run  first  in  a  deterministic 
homogeneous  medium  so  that  results  could  be  easily  interpreted 
and  compared  to  known  methods  and  approximations.  All  test 
cases  are  checked  to  obey  the  criteria  required  by  the  trans¬ 
fer  function. 

The  values  of  the  basic  parameters  such  as  source  depth 
y^,  gradient  g  and  reference  speed  of  sound  c^  were  chosen  in 
order  to  model  the  wave  propagation  between  a  source  located 
on  the  SOFAR  channel  axis  and  a  deep  receiver.  At  moderate 
latitudes  from  60 "S  to  60 the  speed  of  sound  on  the  SOFAR 
axis  ranges  from  1450-1485  m/sec  in  the  Pacific  and  from 
1450-1500  m/sec  in  the  Atlantic  [Ref.  15],  The  depth  of  the 
SOFAR  axis  is  typically  between  1000  and  1200  m  at  moderate 
latitudes  [Ref.  15].  The  gradient  below  the  SOFAR  axis  is 
given  by  g  =  0.017  1/sec  [Ref.  19].  Values  for  the  standard 
deviation  a  of  the  random  component  of  the  index  of  refraction 
range  from  0.25x10  ^  [Refs.  17,18]  to  1.0  x 10  ^  [Ref.  16]. 
These  values  were  used  in  choosing  representative  parameters 
for  the  following  test  cases; 


Homogeneous  Cases 


As  a  baseline  test  case  we  define  the  deterministic 


homogeneous  case  HMGl ; 
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=  1475  m/sec. 

g  =  0:  constant  speed  of  sound, 
a  =  0:  deterministic  medium. 

M  =  11,  N  =  11;  a  11  by  11  element  transmit  array  with 
rectangular  amplitude  weighting  in  both  the  x 
and  y  directions. 

M'=5,N'=5:  a5by5  element  receive  array. 


=  0.0  m 

o 

X  =0.0 
r 

m 

=  1000.0  m 

y  =  2000.0 

m 

o 

r 

=  0.0  m 

o 

z^  =  1732.05 

m 

,  =  d  =  d'  = 

X  y  X 

d'  =  0.7375  m; 

y 

spacing  equal  to 

0  =  30°,  =  90°:  spherical  angles  denoting  the  line 

of  sight  between  the  center  of  the  transmit  array 
and  the  center  of  the  receive  array. 

=  0,  V  =0.5:  line  of  sight  direction  cosine  values 

from  the  center  of  the  transmit  array  to  the  center 
of  the  receive  array. 

%  “  '^b  ”  0*5  5  direction  cosine  values  used  to  steer 

the  main  lobe  of  the  transmit  far-field  beam  pattern 
towards  the  center  of  the  receive  array. 

KMAX  =  1:  single  frequency  component. 

f  =  1000  Hz,  T  =  0.001  sec. 
o  o 

c^  =  aj^  =  1.0;  complex  Fourier  series  coefficient  equals  1. 
L  =  3;  3  times  samples  per  fundamental  period  T^. 

Based  on  the  definition  of  test  case  HMGl  we  define 


the  following 
HMG2 ;  Same 

^r  = 

'^r  = 
9  = 


r 


cases : 

as  HMGl  except 

200  m,  =  1720.464  m 

u,  =  0.1,  V  =v,  =  0.5 

b  r  b 

30.66°,  4/^  =  78.69°. 
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Same  as  HMGl  except 

=  1500  m,  =  866.0246  m 

Same  as  HMGl  except 

=  2500  m,  =  2598.074  m 


Scune 

as  HMGl  except 

^x  = 

d  =  d'  =  d'  = 
y  X  y 

.36875  m 

l-h 

o 

II 

2000  Hz,  T^  = 

=  0 

.0005  sec. 

Same 

as  HMGl  except 

X  = 

r 

1224.745  m 

z 

r 

=  1224.745  m 

Uj^  =  0.6124, 

V 

r 

in 

• 

o 

II 

II 

CD 

II 

52.24° , 

4' 

^r 

=  39.23°. 

Same  as  HMGl  except 


=  d  =  d'  =  d'  =  0.184375  m 
X  y  X  y 


=  4000  Hz, 


Tq  =  0.00025  sec. 


Same  as  HMG2  except 

KMAX  =  5:  five  frequency  components. 


f^  =  1000  Hz, 

=  5000  Hz. 

max 


T  =  0.001  sec. 
o 


fj^  =  3000  Hz:  frequency  used  for  steering  the 
transmit  beam  pattern. 

d  =  d  =  d'  =  d'  =  1475/2  f  =  0.1475  m 
x  y  X  y  max 

c-^  =  a,  =  1.0;  k  =  1,2,...,  5:  all  complex  Fourier 
series  coefficients  equal  1. 

L  =  11:  11  time  samples  per  fundamental  period  T^. 


Deterministic  Inhomocreneous  Cases 


The  following  deterministic  inhomogeneous  cases 


incorporate  a  sound  speed  profile  with  a  constant  gradient. 

The  definitions  are  based  on  those  used  for  the  homogeneous  cases 


INHMGl,  INHMG2,  and  INHMG5 ;  Same  as  HMGl ,  HMG2  and  HMG5  except 
g  =  0.017  1/sec. 

INHMG8 :  Same  as  HMG8  except 

g  =  0.017  1/sec. 

3 .  Random  Inhomogeneous  Cases 

The  following  cases  with  an  index  of  refraction  com¬ 
posed  of  a  deterministic  and  a  random  component  are  defined: 

RINHMGl,  RINHMG2:  Same  as  INHMGl,  INHMG2  except 
a  =  1.0x10  ^ 

RINHMG8:  Same  as  INHMG8  except 

a  =  0.25  X  lO""^ 


C.  RESULTS 

This  section  will  present  and  discuss  the  results  from 
the  computer  model  runs  on  the  test  cases  as  defined  in  Section 
IV. B.  Some  cases  were  run  using  both  methods  of  integration 
w.r.t.  u^.  The  two  methods  showed  no  significant  differences. 

First,  some  trivial  checks  were  made  using  the  various 
homogeneous  cases.  Note,  that  all  test  cases  have  the  same 
line  of  sight  direction  cosine  v^.  For  the  homogeneous  cases 
with  a  single  frequency  component  of  1000  Hz  the  amplitude 
and  phase  terms  of  the  ocean  transfer  function  are  equal  and 
given  by: 


since  those  cases  have  the  same  ocean  transfer  function 


values,  the  difference  in  the  magnitude  of  the  frequency 
response  H(f,i,q)  for  the  various  cases  should  only  depend 
on  the  total  range  given  by 

"total  - 

and  the  transmit  far-field  beam  pattern.  Cases  HMGl,  HMG3 
and  HMG4  also  have  the  same  line  of  sight  direction  cosine 
value  u^,  so  the  influence  of  the  transmit  far-field  beam 
pattern  is  the  scime  for  these  cases  and  magnitude  differences 
between  the  frequency  responses  can  only  depend  on  the  rotal 
range  and  should  behave  in  accordance  with  the  expected 
spherical  spreading  for  a  homogeneous  medium.  Table  IV. 1 
shows  cases  HMGl,  HMG3  and  HMG4  together  with  their  total 
range  and  the  magnitude  of  the  frequency  response  at  the  center 
element  of  the  receive  array,  i.e.,  |H(f, 0,0)1.  We  conclude 
that  the  magnitude  indeed  falls  off  as  1/r,  only  the  total 
range  is  of  importance. 

Table  IV. 2  shows  cases  HMGl,  HMG2  and  HMG6  which  differ 
only  in  cross  range  - x^  and  line  of  sight  direction  cosine 
value  u^  but  have  the  same  total  range.  In  these  cases  only 
the  far-field  beam  pattern  can  cause  differences  in  the  magni¬ 
tude  of  the  frequency  response.  The  teible  shows  that  as  the 
cross  range  increases  the  magnitude  falls  off,  which  is  con¬ 
sistent  with  the  fact  that  when  a  beam  pattern  is  steered 
towards  end-fire  (to  larger  values  of  and  Vj^)  the 


•  .1 
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directivity  index  decreases  due  to  the  increasing  beamwidth 
of  the  main  lobe. 

Table  IV. 3  shows  the  difference  in  magnitude  of  the  fre¬ 
quency  response  at  the  center  of  the  receive  array  between 
HMGl,  HMG5 ,  and  HMG7 .  Cases  with  the  same  geometry  but  a 
different  value  for  the  single  frequency  component.  This 
difference  in  magnitude  is  due  to  the  frequency  dependent 
amplitude  factor  inherent  in  the  WKB  approximation,  which  is 
not  exact  for  a  homogeneous  medium. 

Table  IV. 4  compares  the  magnitudes  of  the  frequency 
response  between  the  homogeneous  cases  HMGl ,  HMG2 ,  and  HMG5 
and  their  inhomogeneous  counterparts  INHMGl,  INHMG2  and  INHMG5. 
We  see  that  the  magnitudes  for  the  inhomogeneous  cases  are 
slightly  lower,  an  expected  phenomena  because  of  the  channel¬ 
ing  effect  in  the  inhomogeneous  medium.  In  a  medium  with  a 
constant  positive  gradient,  the  acoustic  rays  are  bent  upwards 
and  the  soxind  intensity  decreases  with  increasing  depth. 

This  causes  a  lower  magnitude  of  the  acoustic  signal  then 
expected  by  spherical  spreading  only. 

We  will  now  present  numerical  and  graphical  output  from 
the  three  dimensional  DFT  beamforming  program.  The  frequency 
transform  for  the  single  frequency  component  cases  simply 
reveals  the  magnitude  of  the  transfer  function  multiplied  by 
the  amplitude  of  the  input  electrical  signal  and  is  not  shown 
here.  Figures  7  through  10  show  the  magnitude  of  the  spatial 

Ygj^(q,r,n)  versus  direction  cosine  u,  and 


DFT  w.r.t.  X,  i.e.. 


the  magnitude  of  the  spatial  DFT  w.r.t.  y,  i.e.,  Y  ^(q,m,s) 
versus  direction  cosine  u,  for  the  cases  HMGl,  HMG2 ,  INHMGl 
and  INHMG2 .  These  graphs  allow  us  to  estimate  the  direction 
of  propagation  of  the  incoming  plane  wave  by  finding  the 
coordinates  of  the  maxima  on  the  graphs.  Tables  IV. 5  to 
IV. 8  provide  numerical  values  for  the  above  mentioned  DFT ' s 
w.r.t.  X  and  y.  They  show  values  of  Ygjj(q,r,n)  and  Ygjj(q,m,s) 
around  the  maxima  and  make  it  possible  to  form  more  accurate 
estimates . 

From  the  estimated  values  u  and  v  ,  we  can  calculate 

o  o 

estimates  of  target  elevation  angle  0^  and  target  bearing 

angle  Table  IV. 9  gives  the  results  for  the  different 

cases.  The  values  of  the  estimators  u  and  v  were  determined 

o  o 

by  performing  a  second  order  interpolation  on  the  data  given 
in  tables  IV.  5  through  IV. 8.  It  is  seen  that  the  estimates 
of  target  angles  for  the  homogeneous  cases  are  correct  and 
equal  the  line  of  sight  angles  from  the  receiver  to  the  trans¬ 
mitter.  The  estimated  target  angles  for  the  deterministic 
inhomogeneous  cases  deviate  slightly  from  the  line  of  sight 
angles.  This  is  due  to  the  positive  sound  speed  gradient 
which  causes  upward  bending  of  the  acoustic  rays  so  the  direc¬ 
tion  cosine  v^,  directly  related  to  the  vertical  angle  of 
propagation  of  the  plane  wave,  decreases  as  the  wave  travels 
along.  The  net  effect  for  these  cases  is  relatively  small 
but  noticable. 

The  effect  of  randomness,  characterized  by  a  random  component 
in  the  index  of  refraction  with  a  standard  deviation  of 


-4 

a  =  1.0  X 10  is  examined  next.  Figure  11  shows  estimates  of 
u  and  for  the  case  RlNHMGl  and  Fig.  12  shows  estimates  of 
Vq  for  cases  RlNHMGl  and  RINHMG2 .  As  can  be  seen,  the  esti¬ 
mate  of  is  not  influenced,  which  follows  readily  from  the 
fact  that  the  randomness  is  only  present  in  the  y  direction. 

The  estimate  of  as  shown  in  Fig.  12  for  the  case  RINHMG2 
shows  the  same  random  fluctuations  as  the  one  shown  for  case 
RlNHMGl  in  Fig.  11,  since  both  cases  used  the  same  seed  for 
the  random  number  generator,  resulting  in  the  same  random 
number  sequence.  Figure  12  also  shows  RlNHMGl  with  a  different 
seed  for  the  random  number  sequence.  Although  both  are  only 
two  possible  samples  the^  clearly  show  the  strong  effect  of 
the  randomness  on  the  beamformer  output  and  the  deterioration 
of  the  quality  of  the  estimates  (see  Table  IV. 9) . 

Two  three-dimensional  plots  are  shown  in  Fig.  13,  one  for 
the  deterministic  inhomogeneous  case  INHMG2  and  one  for  the 
random  inhomogeneous  case  RINHMG2 .  These  graphs  illustrate 
the  kind  of  output  that  can  be  obtained  from  the  DFT  beam- 
former,  although  it  is  difficult  to  accurately  read  off 
numerical  values. 

Finally,  the  cases  HMG8 ,  INHMG8  and  RINGHM8 ,  which  involve 
several  frequency  components,  are  considered.  The  output 
for  these  cases  is  given  in  numerical  and  graphical  form. 

Table  IV. 10  compares  the  Fourier  series  coefficients  of  the 
output  electrical  signal  with  the  Fourier  series  coefficients 
of  the  transmitted  electrical  signal.  The  widely  varying  values 
are  caused  by: 


1)  The  beam  pattern  of  the  transmit  array  is  frequency 
dependent.  The  differences  between  the  transmitted  frequency 
components  are  relatively  large  and  cause  large  differences  in 
amplitude  c  che  Fourier  series  coefficients.  The  maximum 
values  are  observed  at  f  =  3000  Hz,  i.e.,  the  frequency  for 
(/hich  beam  steering  was  done. 

2)  The  deterministic  inhomogeneous  medium:  As  already 
shown  in  other  test  cases  the  amplitude  of  the  received  signal 
is  frequency  dependent  due  to  the  WKB  approximation  to  the 
wave  equation  in  an  inhomogeneous  medium.  Equation  (2.10) 
clearly  shows  the  frequency  dependence  of  the  medium  transfer 
function . 

3)  The  random  inhomogeneous  behavior  of  the  medium;  The 
randomness  introduces  additional  deviations.  It  can  be  seen 
from  table  IV. 10  that  the  effect  of  the  randomness  increases 
with  increasing  frequency.  Also,  in  the  case  of  a  random 
medium,  the  estimation  of  the  Fourier  series  coefficients 
becomes  strongly  dependent  on  the  location  of  the  element  in 
the  receive  array.  Note  that  table  IV. 10  only  shows  values 
for  the  receive  element  (m  =  0,  n  =  0). 

The  received  angular  spectrum  for  the  different  frequency 
components  of  the  cases  INHMG8  and  RINHMG8  is  shown  in  Figs. 

14  through  19.  For  the  deterministic  inhomogeneous  case. 

Figs.  14  through  16  show  a  regular  behavior  of  the  angular 
spectrum  according  to  the  rectangular  amplitude  weighting 
used  for  the  receive  array.  Estimation  of  u^  is  exact. 


ESTIMRTION  OF  VO 


Fig.  12  Estimation  of  v  for  Cases  RINHMG2  and 
RINHMGl  (Different  Seeds) 


ESTIMATION  OF  UO 


ESTIMATION  OF  VO 


Fig.  9 


Estimation  of  Direction  of 
Propagation  for  Case  INHMGl 


TABLE  IV. 10 


Estimation  of  Fourier  Series  Coefficients 
Ygj^(q,0,0)  for  HMG8,  INHMG8  and  RINHMG8 


Frequency 
(in  Hz) 

HMG8 

INHMG8 

RINHMG8 

1000 

0.0010987 

0.0009705 

0.0010341 

2000 

0.0098126 

0.0098901 

0.0102369 

3000 

0.0209366 

0.0201906 

0.0177316 

4000 

0.0138747 

0.0122052 

0.0087726 

5000 

0.0024539 

0.0033176 

0.0051316 

Note;  The  Fourier  series  coefficients  of  the  trans¬ 
mitted  electrical  signal  are  c,  =  a,  =  1; 
k*l,2,...,5  ^  ^ 


TABLE  IV. 11 

A  /V  A 

Estimated  Values  of  u^,  v  ,  9  and  \b 

o  o  o 


Case  Frequency  u 
(in  Hz)  *■ 


(in  degrees) 


INHMG8 


'Po 

(in  degrees) 


RINHMG8 


5 


TABLE  IV. 9 


Estimated  Values  u  ,  v  ,  0  and  ip 

o  o  o  o 


Case 

^O 

V 

o 

(in  degrees) 

(in  degrees) 

HMGl 

0.000 

0.500 

+30.00 

+90.00 

HMG2 

0.100 

0.500 

+30.65 

+78.69 

INHMGl 

0.000 

0.492 

+29.47 

+90.00 

INHMG2 

0.101 

0.492 

+30.15 

+78.40 

RINHMG2 

0.101 

0.523 

+32.19 

+79.07 

(seed 

#1) 

RINHMGl 

0.000 

0.523 

+31.53 

+90.00 

(seed 

#1) 

RINHMGl 

0.000 

0.508 

+30.53 

+90.00 

( seed 

#2) 

1)  The  line  of  sight  angles  from  receive  to  transmit 
array  (or  target  angles)  are: 

HMGl,  INHMGl,  RINHMGl :  0^.  =  +30.00®,  \p^  =  +90.00® 


HMG2,  INHMG2,  RINHMG2 :  0  =  +30.65°,  ^  =  +78.69® 

r  r 


2)  Note  that  equation  (4.14)  which  is  used  to  calculate 
estimates  of  target  bearing  angle  from  estimates  of 
Uq  and  Vq  returns  two  results.  This  is  due  to  the 
cimbiguity  of  a  planar  array.  For  the  above  table  the 
appropriate  values  were  selected. 


TABLE  IV. 7 


Magnitude  of 
Cosine 

Ysn (q,m,s) 
V  for  HMGl 

Vs.  Direction 
,  HMG2 

Direction 

Cosine  v 

HMGl 

HMG2 

0.484 

.0121367 

.0120534 

0.489 

.0121531 

.0120696 

0.494 

.0121636 

.0120800 

0.499 

.0121681 

.0120844 

0.504 

.0121666 

.0120828 

0.509 

.0121591 

.0120754 

0.514 

.0121457 

.0120620 

TABLE  IV. 8 

Magnitude  of 
Cosine 

YsM(q,in,s)  Vs.  Direction 

V  for  INHMGl,  INHMG2 

Direction 

Cosine  v 

INHMGl 

INHMG2 

0.475 

.0116731 

.0115947 

0.480 

.0116893 

.0116111 

0.485 

.0117001 

.0116215 

0.490 

.0117051 

.0116264 

0.495 

.0117045 

.0116258 

0.500 

.0116983 

.0116196 

0.505 


0116864 


.0126078 


TABLE  IV. 5 


Magnitude  of 
Cosine 

u  for  HMGl, 

Vs.  Direction 
INHMGl 

Direction 

Cosine  u 

HMGl 

INHMGl 

0.015 

.0121416 

.0116804 

0.010 

.0121566 

.0116945 

0.005 

.0121655 

.0117030 

0.0 

.0121685 

.0117058 

0.005 

.0121655 

.0117030 

0.010 

.0121566 

.0116945 

0.015 

.0121416 

.0116804 

TABLE  IV. 6 

Magnitude  of 
Cosine 

YSN(q»r,n) 
u  for  HMG2 , 

Vs.  Direction 
INHMG2 

Direction 

Cosine  u 

HMG2 

INHMG2 

0.085 

.0120593 

.0116000 

0.090 

.0120744 

.0116147 

0.095 

.0120836 

.0116239 

0.100 

.0120869 

.0116274 

0.105 

.0120842 

.0116254 

0.110 

.0120756 

.0116178 

0.115 


0120611 


0116046 


TABLE  IV. 3 


Magnitude  of  Frequency  Response  Comparing  Cases 
Which  Are  Different  Only  in  Frequency 


Case  Frequency  )H(f,0,0)| 


(in  H2) 

HMGl 

1000 

0.02434 

HMG5 

2000 

0.03441 

HMG7 

4000 

0.04866 

TABLE  IV. 4 

Magnitude  of  Frequency  Response  Comparing 
Homogeneous  and  Inhomogeneous  Cases 


Case 

H(f  ,0 

,0) 

HMGl  - 

INHMGl 

0.02434 

-  0.02341 

HMG2  - 

INHMG2 

0.02417 

-  0.02326 

HMG5  - 

INHMG5 

0.03441 

-  0.03310 

TABLE  IV. 1 


Magnitude  of  Frequency  Response  Comparing  Cases 
Which  Are  Different  Only  in  Total  Range 


Case 

Total  Range 
(in  m) 

|H(f ,0,0) 

HMGl 

2000 

0.02434 

HMG3 

1000 

0.04868 

HMG4 

3000 

0.01622 

TABLE  IV. 2 

Magnitude  of  Frequency  Response  Comparing  Cases 
Which  Are  Different  Only  in  Cross  Range 


Case 


Total  Range 


Cross  Range 


H(f ,0,0) 


0.0 

200.0 


0.02434 

0.02417 


1224.75 


0.01654 


Estimation  of  returns  a  slightly  different  value  from  the 

line  of  sight  value  for  direction  cosine  v^  due  to  the 

inhomogeneous  behavior  of  the  medium.  Figures  14,  15  and 

16  demonstrate  the  increasing  directivity  of  the  receive 

array  with  increasing  frequency,  thus  leading  to  more  accurate 

estimates  for  higher  frequencies. 

Figures  17,  18  and  19  clearly  show  the  strong  influence 

of  the  randomness  on  the  resulting  angular  spectrum.  The 

randomness  creates  strong  "false"  sidelobes  and  causes  an 

offset  of  the  estimated  value  of  direction  cosine  v  .  The 

o 

effects  of  the  randomness  increase  with  increasing  frequency. 

At  f  =  5000  Hz  (Fig.  19)  the  estimate  of  v^  is  determined  by 

a  spurious  side  lobe  and  returns  a  totally  wrong  value. 

Table  IV. 11  shows  estimated  values  of  u  and  v  at  the  differ- 

o  o 

ent  frequencies  and  their  associated  estimates  of  target 

/s 

elevation  and  bearing  angles,  9  and  tp  ,  respectively. 


Fig.  14  Angular  Spectrum  for  Case  INHMG8 
at  f  =  1000  and  f  -  2000  Hz 


Computer  simulation  of  a  mathematical  model  of  wave  propa¬ 
gation  through  a  random,  space-variant,  time-invariant,  ocean 
medium  between  two  planar  arrays  has  been  accomplished.  The 
computer  simulation  incorporates  an  index  of  refraction  which 
is  only  a  function  of  depth.  The  deterministic  component  of 
the  index  of  refraction  represents  a  sound  speed  profile  with 
constant  gradient.  The  random  component  is  assumed  to  have 
statistical  properties  independent  of  depth.  The  effects  of 
the  randomness  on  the  electrical  output  signals  at  the  differ¬ 
ent  receiver  elements  are  assumed  to  be  uncorrelated. 

The  computer  program  was  written  in  FORTRAN.  It  implements 
both  analytical  and  numerical  integration  techniques.  The 
results  from  both  techniques  are  in  close  agreement. 

A  number  of  test  cases  have  been  run.  Interpretation  of 
the  data  for  simple  cases  show  results  which  are  consistent 
with  expectations  from  theory.  The  model  correctly  demonstrated 
the  influences  of  the  far  field  beam  pattern  of  the  transmit 
array,  the  source-receiver  geometry,  and  the  deterministic 
inhomogeneous  behavior  of  the  medium.  Also,  a  random  inhomo¬ 
geneous  medium  was  shown  to  affect  the  wave  propagation  and 
distort  the  frequency  and  angular  spectra,  leading  to  errors 
in  the  estimation  of  the  Fourier  series  coefficients,  target 
bearing  angle,  and  target  elevation  angle.  A  final  test  case 
RINHMG8,  with  a  transmit  signal  containing  several  frequency 


ccmponents  demonstrated  a  complex  dependence  on  the  above 
mentioned  factors,  making  interpretation  of  the  results  diffi¬ 
cult.  The  effect  of  the  randomness  was  very  pronounced  in 
this  case. 

The  computer  program  could  be  further  improved  by  imple¬ 
menting  a  more  sophisticated  method  to  simulate  the  randomness 
of  the  medium.  The  assumption  that  the  random  phase  terms  at 
the  different  receiver  elements  are  uncorrelated  breaks  down 
for  small  interelement  spacings.  The  computer  program  can 
also  be  easily  modified  to  accommodate  a  more  complex  sound 
speed  profile  in  which  the  gradient  is  a  function  of  depth, 
for  example. 

Future  application  of  the  model  could  be  devoted  to 
investigating  the  propagation  of  more  complex  transmit  signals 
One  could  also  investigate  the  effects  of  the  randomness  of 
the  ocean  medium  on  signal  processing  using  different  sizes 
of  receive  arrays,  i.e.,  increasing  the  number  of  elements. 


LIST  OF  REFERENCES 


1.  Ziomek,  L.J.,  Underwater  Acoustics — A  Linear  System  Theory 
Approach,  Academic  Press,  Orlando,  Florida  (unpublished) . 

2.  Ibid.,  Chapter  7,  Section  7.2.1. 

3.  Ibid.,  Chapter  7,  Section  7.2.2. 

4.  Ibid.,  Chapter  5,  Section  5.1. 

5.  Laval,  R.  ,  ‘'Sound  Propagation  Effects  on  Signal  Processing, 
Signal  Processing,  edited  by  Griffiths,  P.L.  and  others, 
pp.  223-241,  Academic  Press,  New  York,  1973. 

6.  Laval,  R. ,  "Time-Frequency-Space  Generalized  Coherence 
and  Scattering  Functions , "  Aspects  of  Signal  Processing, 
edited  by  Tacconi ,  G. ,  Vol .  I,  pp.  69-87,  D.  Reidel 
Publishing  Company,  Dordrecht,  The  Netherlands,  1977. 

7.  Laval,  R.  and  Labasque ,  Y.,  "Medium  Inhomogeneities  and 
Instabilities;  Effects  on  Spatial  and  Temporal  Processing, 
Underwater  Acoustics  and  Signal  Processing,  edited  by 
Bjorno,  L. ,  pp.  41-70,  D.  Reidel  Publishing  Company , 
Dordrecht,  The  Netherlands,  1981. 

8.  Middleton,  D.,  "A  Statistical  Theory  of  Reverberation  and 
Similar  First-Order  Scattered  Fields.  Part  III:  Waveforms 
and  Fields,"  IEEE  Trans.  Information  Theory,  Vol.  18, 

pp.  35-67,  1972. 

9.  Middleton,  D.,  "The  Underwater  Medium  as  a  Generalized 
Communication  Channel,"  Underwater  Acoustics  and  Signal 
Processing ,  edited  by  Bjorno,  L. ,  pp.  589-612,  D.  Reidel 
Publishing  Company,  Dordrecht,  The  Netherlands,  1981. 

0.  Officer,  C.B.,  Introduction  to  the  Theory  of  Sound  Trans¬ 
mission,  pp.  67-68,  McGraw-Hill,  New  York,  1958. 

1.  Gerald,  C.T.,  Applied  Numerical  Analysis,  pp.  60-72, 

Addison -Wes ley ,  1970  . 

2.  ^Moursund,  D.G.  and  Duris,  C.S.,  Elementary  Theory  and 
Application  of  Numerical  Analysis,  pp.  194-199,  McGraw- 
Hill,  1967. 

3.  Squire,  W. ,  Integration  for  Engineers  and  Scientists, 
American  Elsevier  Publishing  Company,  New  York,  1970. 


Brekhovskikh ,  L.M.,  Waves  in  Layered  Media,  pp.  234-239, 
Academic  Press,  New  York,  1980. 

Brekhovskikh,  L.  and  Lysanov,  Yu.,  Fundamentals  of  Ocean 
Acoustics ,  p.  3,  Springer-Verlag ,  1982. 

Clarke,  R.H.,  "Sound  Propagation  in  a  Variable  Ocean," 
Journal  of  Sound  and  Vibration,  Vol.  34,  p.  472,  1974. 

Mintzer,  D.,  "Wave  Propagation  in  a  Randomly  Inhomogeneous 
Medium  I,"  Journal  of  the  Acoustic  Society  of  America, 

Vol.  25,  p.  925,  1953. 

Mintzer,  D. ,  "Wave  Propagation  in  a  Randomly  Inhomogeneous 
Medium  I  and  II,"  Journal  of  the  Acoustic  Society  of 
America ,  Vol.  26,  p.  1110,  1953. 

Kinsler,  L.E.,  and  others.  Fundamentals  of  Acoustics , 
p.  401,  3rd  ed.,  Wiley  &  Sons,  1981. 


INITIAL  DISTRIBUTION  LIST 


No 


1.  Defense  Technical  Information  Center 
Cameron  Station 

Alexandria,  Virginia  22314 

2.  Library,  Code  0142 
Naval  Postgraduate  School 
Monterey,  California  93943 

3.  Department  Library,  Code  61 
Department  of  Physics 
Naval  Postgraduate  School 
Monterey,  California  93943 

4.  Assistant  Professor  L.J.  Ziomek,  Code  62Zm 
Department  of  Electrical  and 

Computer  Engineering 
Naval  Postgraduate  School 
Monterey,  California  93943 

5.  Inspecteur  Onderwijs  Zeemacht 
Ministry  of  Defense  (Navy) 

Koningin  Marialaan  17 

Den  Haag 

The  Netherlands 

6.  Vlagofficier  belast  met  de  of ficiersvorming 
Koninklijk  instituut  voor  de  Marine 

Het  Nieuwe  Diep  8 

Den  Helder 

The  Netherlands 

7.  LT  J.  Vos 

p/a  Dorpsstraat  12 
Schipluiden,  2636  CC 
The  Netherlands 

8.  Mr.  Charles  Stewart 
DARPA 

1400  Wilson  Blvd. 

Arlington,  VA  22200 


Copies 

2 

2 

1 

15 

2 


2 


3 


1 


FILMED 

6-85 


DTIC 


